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Abstract 

In the multidimensional CESE development, triangles and tetrahedra turn out to be the most natural 
building blocks for 2D and 3D spatial meshes. As such the CESE method is compatible with the simplest 
unstructured meshes and thus can be easily applied to solve problems with complex geometries. However, 
because the method uses space-time staggered stencils, solution decoupling may become a real nuisance in 
applications involving unstructured meshes. In this paper we will describe a simple and general remedy 
which, according to numerical experiments, has removed any possibility of solution decoupling. Moreover, in 
a real-world viscous flow simulation near a solid wall, one often encounters a case where a boundary with high 
curvature or sharp corner is surrounded by triangular/tetrahedral meshes of extremely high aspect ratio (up 
to 10 6 ). For such an extreme case, the spatial projection of a space-time compounded conservation element 
constructed using the original CESE design may become highly concave and thus its centroid (referred to as a 
spatial solution point) may lie far outside of the spatial projection. It could even be embedded beyond a solid 
wall boundary and causes serious numerical difficulties. In this paper we will also present a new procedure 
for constructing conservation elements and solution elements which effectively overcomes the difficulties 
associated with the original design. Another difficulty issue which was addressed more recently is the well- 
known fact that accuracy of gradient computations involving triangular/tetrahedral grids deteriorates rapidly 
as the aspect ratio of grid cells increases. The root cause of this difficulty was clearly identified and several 
remedies to overcome it were found through a rigorous mathematical analysis. However, because of the 
length of the current paper and the complexity of mathematics involved, this new work will be presented in 
another paper. 


1. Introduction 

The space-time conservation element and solution element (CESE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-70]. Its nontraditional features include: 
(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a time marching strategy that 
has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated without using 
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any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability to capture shocks 
without using Riemann solvers, and also partially contributes to the unusually simple CESE non-reflecting 
boundary conditions); (iv) the requirement that each scheme be built from a nondissipative core scheme and, 
as a result, the numerical dissipation can be controlled effectively; and (v) the fact that mesh values of the 
physical dependent variables and their spatial derivatives are considered as independent mesh variables to be 
solve for simultaneously. Note that CEs are non-overlapping space-time subdomains introduced such that (i) 
the computational domain can be filled by these subdomains; and (ii) flux conservation can be enforced over 
each of them and also over the union of any combination of them. On the other hand, SEs are space-time 
subdomains introduced such that (i) the boundary of each CE can be divided into component parts with 
each of them belonging to a unique SE; and (ii) within a SE, any physical flux vector is approximated using 
simple smooth functions constructed from the solution information stored at a space-time solution point. In 
general, a CE does not coincide with a SE. 

Without using flux-splitting, dimensional-splitting, mesh-alignment or other special techniques, since its 
inception in 1991 [1], the unstructured-mesh compatible 2nd-order CESE method has been used to obtain 
numerous accurate ID, 2D and 3D steady and unsteady flow solutions with Mach numbers ranging from 
0.0028 to 10 [15]. The physical phenomena modeled include traveling and interacting shocks, acoustic waves, 
vortex shedding, viscous flows, detonation waves, cavitation, flows in fluid film bearings, heat conduction 
with melting and/or freezing, MHD vortex, hydraulic jump, crystal growth, chromatographic problems, 
solar wind and stress waves in elastic solids [3-70]. In particular, its unexpected simple non-reflecting 
boundary conditions [12,50] and rather unique capability to resolve both strong shocks and small disturbances 
(e.g., acoustic waves) simultaneously [19,21,22,66,67] make the CESE method an effective tool for attacking 
computational aeroacoustics (CAA) problems. Note that the fact that the second-order CESE method is 
capable of solving CAA problems accurately is an exception to the commonly-held belief that a second-order 
scheme is not adequate for modeling CAA problems. Also note that, while numerical dissipation is needed 
for shock capturing, it may also result in annihilation of small disturbances. Thus a solver that can handle 
both strong shocks and small disturbances simultaneously must be able to overcome this difficulty. 

In spite of its nontraditional features and robust capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the results of a pursuit driven by these simple ideas. The first 
and foremost is the belief that the method must be solid in physics. As such, in the CESE development, 
conservation laws are enforced locally and globally in their natural space-time unity forms for ID, 2D and 3D 
cases. Moreover, because direct physical interaction generally occurs only among the immediate neighbors, 
use of the simplest stencil also becomes a common CESE feature. Obviously, this requirement has the effect 
of simplifying boundary-condition implementation. 

The second idea emerges from the realization that stability and accuracy are two competing issues in 
time-accurate computations, i.e., too much numerical dissipation will degrade accuracy while too little of 
it will cause instability. In other words, to meet both accuracy and stability requirements, computation 
must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a 
real dilemma in numerical method development. As an example, high order schemes generally have higher 
accuracy and lower numerical dissipation. However, they are susceptible to computational instabilities. In 
fact, in complicated real-world applications, not only they seldom live up to their nominal order of accuracy — 
generally they possess only first-order accuracy when shocks are present, stability of these schemes often is 
highly problem-dependent and difficulty to predict. Moreover, generally speaking, the higher the order of 
the scheme, the more restritive its stability conditions become. To confront this issue head-on, in the CESE 
development, generally it is required that a solver be built from a nondissipative (i.e., neutrally stable) core 
scheme. By definition, computations involving a neutrally stable scheme are performed right on the edge of 
instability and therefore the numerical results generated are nondissipative. As such numerical dissipation 
can be controlled effectively if the deviation of a solver from its nondissipative core scheme can be adjusted 
using some built-in parameters. 

Moreover, because an accurate viscous flow simulation requires that the numerical dissipation be much 
smaller than the physical dissipation which decreases as the Reynolds number increases, in principle, an 
accurate and robust solver for high-Reynolds-number flows must be able to cut numerical dissipation as the 
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Reynolds number increases. Obviously this requirement can only be met by a solver built from a 
nondissipative core solver. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 

(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of ID 
schemes so that multidimensional shocks can be captured without using any mesh alignment technique ; 

(iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the method is 
compatible with the simplest unstructured meshes and thus can be easily used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, realization of the above lesser ideas (i)-(iv) follows naturally from 
the first two core ideas. 

Building on its initial successes, efforts to refine and improve the CESE method have continued in the 
past few years. As an example, it was shown in [5] that the numerical dissipation of a dissipative extension 
of a CESE core scheme may increase to an intolerable level as the value of the CFL number decreases from 
its maximum stability bound. As such, in a case with a large CFL number disparity (e.g., a simulation 
with a highly non-uniform spatial mesh and a spatially independent time step), the performance sensitivity 
with respect to the CFL number can lead to a solution that is highly dissipative in a region where the 
local CFL number <C 1. Even though a remedy was suggested in [5], a simple and robust solution to this 
problem had not arrived until a family of new Courant number insensitive schemes [25,31,46-49,51,54,55,58] 
was developed with fresh insights. Note that these new schemes have one important advantage, i.e., all 
variables at each mesh point can be evaluated explicitly without resorting to solving a system 
of linear/nonlinear algebraic equations involving these local mesh variables even in applications 
where systems of multidimensional nonlinear PDEs are solved. 

As another example, even though they are accurate enough to solve CAA problems, second-order CESE 
solvers are not capable of resolving fine flow structures within a boundary layer without using relatively fine 
meshes and/or meshes with large aspect ratios. To overcome this limitation, two neutrally stable 4th-order 
schemes, referred to as the a(3) and a-4 schemes [59-61], were developed in 2006 and intended to serve 
as the core schemes of other high order CESE schemes. Unfortunately, the dissipative extensions of these 
new schemes turned out to have the same intractable stability problem that afflicts traditional high order 
schemes. With the aid of a conceptual leap, this difficulty was finally overcome through the development 
of a new approach [63] by which one can build from a given 2nd-order CESE scheme its 4th-, 6th-, 8th-,... 
order versions which have the same stencil and same stability conditions of the 2nd-order scheme, 
and also retain all other advantages of the latter scheme. As such, these new high order schemes 
can avoid the common shortcomings of traditional high order schemes including: (a) susceptibility to 
computational instabilities; (b) computational inefficiency due to their local implicit nature 
(i.e., at each mesh points, need to solve a system of linear /nonlinear equations involving all the 
mesh variables associated with this mesh point); (c) use of large and elaborate stencils which 
complicates boundary treatments and also makes efficient parallel computing much harder; 
(d) difficulties in applications involving complex geometries; and (e) use of problem-specific 
techniques which are needed to overcome stability problems but often cause undesirable side 
effects. 

In the current paper, we will describe other recent developments which successfully address some of 
the few remaining CESE issues. Some of the new techniques devised have been implemented in the NASA 
code ez4d (developed by Dr. Chau-Lyan Chang, the second author) and the Jacobs Technology Inc. code 
JUSTUS (developed by Dr. Joseph C. Yen, the third author). Their effectiveness has been verified by 
applications involving hypersonic viscous flow over rough solid surfaces [62,64,65] and also by applications 
involving studies of the WICS (weapons Internal Carriage and separation) test data produced at Arnold 
Engineering Development Center [66,67]. Note that the current paper would represent the first systematic 
documentation of the new techniques even though the improved numerical results were presented in [62,64- 
67]. To explain clearly these new techniques and their significance, we will begin our discussion by providing 
some basic information about the CESE method. 
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For simplicity, first we review existing 2nd-order CESE solvers for the simple partial differential equation 
(PDE) 

dii du 

a— = 0 (1.1) 


dt 


dx 


where the advection speed a ^ 0 is a constant. Let x\ = x, and X 2 = t be the coordinates of a two-dimensional 
Euclidean space E 2 . Then, because Eq. (1.1) can be expressed as V • h = 0 with h = (au,u) (i.e., au and u 
are the X\- and x 2 - components of h, respectively), Gauss’ divergence theorem in the space-time E 2 implies 
that Eq. (1.1) is the differential form of the integral conservation law 


(f h ■ ds = 0 (1.2) 

Js(v) 

As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time region V in E 2 , and (ii) 
ds = da n with da and n, respectively, being the area and the unit outward normal vector of a surface 
element on S(V). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through 
the surface element ds, Eq. (1.2) simply states that the total space-time flux of h leaving V through S(V) 
vanishes; (ii) in E 2 , da is the length of a differential line segment on the simple closed curve S(V); and (iii) 
all mathematical operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean 
space. 

To proceed, let Q -| denote the set of all space-time staggered mesh points (j. n) in E 2 (dots in Fig. 2(a)), 

where n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n ± 1/2, n ± 3/2, n ± 5/2, Thus (j ± n) is an 

half integer for each (j, n) € fii. Each (j, n) G fii is associated with a solution element, i.e., SE (j,n). By 
definition, SE(j, n) is the interior of the space-time region bounded by a dashed curve depicted in Fig. 2(b). It 
includes a horizontal line segment, a vertical line segment, and their immediate neighborhood. By definition, 
the end points of the line segments referred to above are excluded from SE(j, n) so that two SEs will not 
overlap. 

Eq. (1.2) will be simulated numerically assuming that, for any ( x,t ) € SE(j, n), u{x,t ) and h(x,t), 
respectively, are approximated by 

u*(x,t ; j,n) = f u] ± (u x )](x - Xj) + {u t )](t - t n ) (j,n) G fii (1.3) 

and 

h*(x,t;j,n) = f ( au*(x,t;j,n ), u* (x,t ; j,n)) ( x,t ) £ SE (j,n) and ( j,n ) G Oi (1.4) 

Note that: (i) u", (iix)™, and (m*)" are constants in SE(j, n), and (in a rough sense) they can be considered 
as the numerical analogues of the values of it, du/ dx, and du/dt at the mesh point ( j,n ), respectively, (ii) 
( Xj,t n ) are the coordinates of the mesh point ( j,n ) with Xj = j Ax and t n = nAt, and (iii) Eq. (1.4) is the 
numerical analogue of the definition h = (ait, it). 

Let u = u*(x,t;j,n) satisfy Eq. (1.1) within SE(j, n). Then one has 

(ut)j = ~a(u x )j (, j,n ) G fii (1.5) 


As a result, Eq. (1.3) reduces to 

u*(x, t ; j, n) = u" ± (u x )” [(x — Xj) — a(t — t n )] (j,n) G fii (1.6) 

i.e., Uj and (it x )" are the only independent mesh variables associated with ( j,n ). 

Let E 2 be divided into non-overlapping rectangular regions (see Fig. 2(a)) referred to as basic conser- 
vation elements (BCEs). As depicted in Figs. 2(c)-2(e), (i) each (j, n) G Oi is assigned with two BCEs, 
i.e., CE_(j, n) and CE + (j, it); (ii) each BCE has one and only one pair of diagonally opposite vertices which 
belong to ; (iii) the space-time E 2 can be filled by CE_(j, n) and CE + (j, n), (j, n) G fii; and (iv) CE(j, n), 
which is the union of CE_(j, n) and CE + (j, n), is referred to as a compounded conservation element (CCE). 
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Note that, among the line segments forming the boundary of CE_ (j, n), AB and AD belong to SE(j, n), 
while CB and CD belong to SE(j — 1/2, n— 1/2). Similarly, among the line segments forming the boundary 
of CE + (j, n), AF and AD belong to SE (j, n), while EF and ED belong to SE(j + 1/2, n — 1/2). As will be 
shown, by imposing the two local flux conservation conditions at each (j, n) G Hi, i.e., 


o 

II 

Tco 

* 

^0- 

(j,n) G Hi 

(1.7) 

/S(CE + (j,n)) 

o 

II 

tec 

* 

^e- 

(j, n) G Hi 

(1.8) 

/s(CE_o») 


one can obtain two equations for two unknowns u ” and {u x )^. 

Let 

v = a — and (u $ ) = — {u x ) j (j, n) G Hi (1.9) 

Ax J 4 J 

where v is the Courant number, and (u®)" is the normalized version of (uj,)". Then, with the aid of Eqs. (1.4), 
(1.6) and (1.9), it can be shown that Eqs. (1.7) and (1.8) are equivalent to 

(1 - v) [u + (1 + v)u x ]™ = (1 - v) [u - (1 + ^)ux]"+i/ 2 2 (j, n) G Hi (1-10) 

and 

(1 + v) [u - (1 - y)%]" = (1 + u) [u + (1 - (j,n)G Oi (1.11) 

respectively. To simplify notation, in the above and hereafter we adopt a convention that can be explained 
using an expression on the left side of Eq. (1.10) as an example, i.e., 

[u + (1 + v)u x ]™ = u] + (1 + v)(u x )j 

By adding Eqs. (1.10) and (1.11) together, one has 

u] = l, {(1 ~ v) [u - (1 + + (1 + u) [u + (1 - } {j, n) G fb (1.12) 

Let 1 — v 2 0, i.e., 1 — v ^ 0 and 1 + v ^ 0. Then Eqs. (1.10) and (1.11) can be divided by (1 — v) and 
(1 + v), respectively. By subtracting the resulting equations from each other, one has 

(Ux)j = \ { [u - (1 + iy)u x ]^~l - [u + (1 - iz)u x ]^Zl /2 } 0/ n) G Oi (1.13) 

Hereafter, by definition the a scheme is formed by Eq. (1.12) and (1.13) for any v. According the above 
derivation, a solution to the a scheme will satisfy Eqs. (1.10) and (1.11) and therefore the two conservation 
conditions Eqs. (1.7) and (1.8) for any v. However, because Eqs. (1.7) and (1.8) => (i.e., “imply”) Eq. (1.12) 
for any u, but they =£> Eq. (1.13) only if u 2 ^ one concludes that a solution to Eqs. (1.7) and (1.8) may 
not be a solution to the a scheme if v 1 = 1. As such the a scheme represents a stronger condition than that 
of Eq. (1.7) and (1.8). 

The a scheme is an explicit and neutrally stable (i.e., non-dissipative) solver of Eq. (1.1) in its stability 
domain \v\ < 1 [63]. Also, according to Eqs. (1.12) and (1.13), it has a space-time staggered stencil, i.e., 
any (j,n) G Hi is associated with a stencil formed by the three mesh points (j, n), ( j — 1/2, n — 1/2) and 
{j + 1/2, n — 1/2). 

Eqs. (1.7) and (1.8), i.e., the flux conservation conditions over the two BCEs CE_(j, n) and CE + (j, n) 
for all (j, n) G Hi, are enforced by any solution to the a scheme. Because (i) the vector h* at any surface 
element lying on any interface separating two neighboring BCE is evaluated using the information from a 
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single SE, and (ii) the unit outward normal vector on the surface element pointing outward from one of these 
two neighboring BCEs is exactly the negative of that pointing outward from another BCE, one concludes that 
the flux leaving one of these BCEs through the interface is the negative of that leaving another BCE through 
the same interface. As a result of this flux cancellation, the local flux conservation relations Eqs. (1.7) and 
(1.8) lead to a global flux conservation relation, i.e. , the total flux of h* leaving the boundary of any space- 
time region that is the union of any combination of BCEs will also vanish. In particular, because CE (j, n) 
is the union of CE_(j, n) and CE + (j, n), 


® h* ■ ds = 0 (j,n) £ Oi (1.14) 

Js{ CE(,») 

must follow from Eqs. (1.7) and(1.8). In fact, it can be shown that Eq. (1.14) is equivalent to Eq. (1.12). 

There is a family of the dissipative extensions of the a scheme in which only the less stringent conservation 
condition Eq. (1.14) is assumed [5,26,31,63]. Because Eq. (1.14) is equivalent to Eq. (1.12), for each of these 
extensions, it” is still evaluated using Eq. (1.12) while (u s )” is evaluated in terms of the mesh variables at 
( j — 1/2, n — 1/2) and ( j + 1/2, n — 1/2) using an equation different from Eq. (1.13). In fact, (i) the a-e 
scheme is formed by Eq. (1.12) and 

(%)" = \ {[(! - e)u + (2e - 1 - v) «x]"+i/ 2 2 ~ K 1 - e)u + (1 - v - 2 } (j, n) £ fii (1.15) 

where e is a real parameter, and (ii) the c-r scheme is formed by Eq. (1.12) and 

( u x)j = ■ 2 (1 1 +t) { l u ~ ( 1 + 2v ~ r ) u *]"+ 1/2 - [u + (1 - 2v - t)u^Zi /2 } (i. n ) e °i (i- 16 ) 

where r is a real parameter ^ — 1. 

Note that: 

(a) By comparing Eqs. (1.13) and (1.15), it is seen that the a scheme is the special case of the a-e scheme 
when e = 0. 

(b) By using Jordan canonical form theorem [71], it can be shown that the a-e scheme is von Neumann 
stable <t=> (i.e., “if and only if”) either 

0 < e < 1 and |^| < 1 (1.17) 


or 

e=M = l (1-18) 

Also one can show that the a-e scheme becomes more dissipative as the value of e increases from 0 to 1 

[5]- 

(c) By using Jordan canonical form theorem [71], it is shown rigorously in [46] that the c-r scheme is von 
Neumann stable <t=> 

v 2 < 1, t > t 0 (v 2 ), and (^ 2 ,r) ^ (1,1) (1.19) 


where 



'0 

4-s-2,/2(2-s-s 2 ) 

s 

s-l + Vl-2s + 5s 2 


K 2s 

It can be shown that: (i) r 0 (s) is continuous at s = 0, i.e., 

lim t 0 (s ) = t o (0) = 


if s = 0 


if 0<s<A- 
if n < s < 1 


0 


(1.20) 


( 1 . 21 ) 
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(ii) r 0 (s) is consistently defined at s = 3/11, i.e. , both the second and third expressions on the right 
side of Eq. (1.20) implies that r 0 (3/ll) = 1/3; (iii) t 0 (s) is continuously differentiable at s = 3/11, i.e., 


lim r/(s) = lim + r'(s) = 121/90 


(1.22) 


where r'(s) = dT 0 ( s )/ ds \ (iv) t 0 (s) is strictly monotonically increasing in the interval 0 < s < 1; (v) 
r 0 ( 1) = 1; and (vi) 

s < t 0 (s) < , 0 < s < 1 (1-23) 

Moreover, (i) the stability conditions Eq. (1.19) agree completely with those generated through numerical 
experiments [31] (Note: the reader is warned that the condition t > T 0 (y 2 ) that appears in Eq. (1.19) 
here is expressed as r > r 0 (|^|) in Eq. (3.12) of [31], i.e., the function r Q defined here is different from 
the function r 0 introduced in [31]), and (ii) for any given fixed value of \v\ < 1, the c-r scheme tends to 
become more dissipative as the value of r increases from its minimum stability bound t 0 (v 2 ). 

The a scheme and its extensions were defined over fii, the set of space-time staggered mesh points marked 
by dots in Fig. 2(a). Obviously, the a scheme can also be defined over fl 2 , the set of all space-time staggered 
mesh points ( j , n) in E 2 where n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n, n ± 1, n ± 2, n ± 3, — 
By definition, each (j, n) G 0 2 is associated with an whole-integer value of ( j + n), and represented by a 
mesh point not marked by dots in Fig. 2(a). Because one can carry out the time marching (using the a 
scheme or any one of its ID extensions) exclusively over either Oi or fl 2 , there is no need to carry out a 
marching over both Hi and fi 2 . Note that this conclusion generally remains valid even if an ID nonuniform 
spatial mesh is used [31]. The exception may occur for a marching scheme where the time-step size varies 
from one spatial region to another such as the local time stepping procedure described in [42] . 

Note that the two disjoint sets Hi and fl 2 have the property that, if a space-time mesh point (j. n) 
belongs to fii (H 2 ), then its four immediate space-time neighbors (j ±1/2, n) and (j, n± 1/2) must belong to 
n 2 (fii). As a result, for a scheme with a space-time staggered stencil such as the a scheme or the classical 
leapfrog scheme, this implies that a time marching carried out over Hi is completely decoupled from that 
carried out over fl 2 . As such one only needs to carry out time marching over one of these two sets. For 
this reason, the problem of “solution decoupling” (which is famously associated with the leapfrog scheme) 
usually can be avoided in the ID case. 

However, as will be shown shortly, for a space-time mesh built from a multidimensional unstructured 
spatial mesh, one generally cannot define two disjoint sets of space-time mesh points which possess the same 
property of fii and H 2 discussed above. As such, the problem of solution decoupling may become a real 
nuisance in multidimensional CESE simulations involving unstructured spatial meshes. 

Is was shown earlier that, for an ID CESE scheme, (i) each mesh point is assigned with two BCEs; and 
(ii) for each dependent physical variable u, each mesh point (j, n ) is assigned with two independent mesh 
variables u™ and (ua,)™. Similarly [12,13], for a 2D (3D) CESE scheme, (i) each space-time solution point 
(to be defined shortly) is assigned with three (four) BCEs and one CCE, with the CCE being the union 
of the three (four) BCEs; and (ii) for each dependent physical variable, each space-time solution point is 
assigned with three (four) independent mesh variables — note that there are two (three) independent spatial 
derivatives in a problem with two (three) spatial dimension. Because a triangle and a tetrahedron have three 
and four sides, respectively triangles and tetrahedra, respectively, become the simplest and most natural 
building blocks of the spatial meshes for 2D and 3D CESE schemes [10-13,16]. As such, in the following 
discussion about solution decoupling in a CESE simulation involving an unstructured 2D spatial domain, we 
consider an unstructured domain built from triangles. 

Consider a spatial domain formed by triangles of arbitrary shapes (see Fig. 3). Here (i) any two 
neighboring triangles share a common side, and (ii) the vertices and centroids of triangles are marked by 
dots and circles, respectively. As shown in Fig. 3, the triangle A FBD has three neighboring triangles. Let 
(i) point G be the centroid of A FBD, and (ii) points A, C and E, respectively, be the centroids of the three 
neighbors of A FBD. Then points A, B , C, D , E, and F form a hexagon. The centroid of the hexagon 
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is denoted by G' and marked by a cross. Hereafter point G' is referred to as the spatial solution point of 
A FBD. In a similar manner, we define spatial solution points A' , C' , E ' ,... of other triangles. 

A space-time solution point is a mesh point which resides at a time level n = 0, 1/2, 1, 3/2, 2, ... (t = 
t n = nAt at the ?rth time level) and has a spatial solution point being its spatial projection. Given any 
spatial solution point G' and any time level n, let (G',n) denote the space-time solution point which resides 
at the nth time level and has point G' being its spatial projection. As explained in [13], the three BCEs 
associated with ( G',n ) are constructed such that their spatial projections are the quadrilaterals GFAB, 
GBCD , and GDEF, respectively. Also, for a reason explained in [13], independent mesh variables are 
stored at space-time solution points. 

As shown in [13], in a 2D CESE scheme using a triangular spatial mesh, the mesh variables at (G', n ) are 
evaluated in terms of those at ( A n — 1/2), (C , n — 1/2) and (E ' , n — 1/2), i.e., the scheme has a space-time 
staggered stencil formed by one point at the nth time level and three points at the (n — l/2)th time level. 
Also it was shown in [13] that, given an unstructured triangular spatial mesh, the time marching cannot be 
carried out independently over two disjoint sets of space-time solution points unless the triangles forming 
the spatial mesh can be divided into two disjoint sets ’hi and d '2 such that two triangles sharing a common 
side always belong to two different sets. As will be shown, because point F is the common vertex of an odd 
number (5) of triangles, the triangles depicted in Fig. 3 cannot be divided into 'h i and T 2 . 

To prove the last proposition by contradiction, let A FBD £ T i . Then because A FHB and A FBD 
share a common side, by definition, A FHB £ T 2 . Similarly, because (i) AFNH and A FHB share a 
common side, (ii) A FMN and AFNH share a common side, and (iii) A FDM and A FMN share a 
common side, one concludes that AFNH £ $ 1 , AFMN £ '[' 2 , and A FDM £ T 1 . As such, we have arrive 
at the contradiction that both A FBD and A FDM belong to Ti in spite of the fact that they share a 
common side. Because we will reach a similar contradiction by assuming A FBD £ T 2 , the proposition has 
been proved. 

The triangles depicted in Fig. 3 cannot be divided into two disjoint sets and As such, according 
to another proposition stated earlier, the CESE time marching cannot be carried out independently over two 
disjoint sets of space-time solution points constructed from the triangular spatial mesh depicted in Fig. 3. 
Moreover, because the incidence that an odd number of triangles sharing a common vertex occurs regularly 
in an unstructured triangular mesh, the CESE time marching generally cannot be carried out independently 
over two disjoint sets of space-time solution points if they are constructed from an unstructured triangular 
spatial mesh. 

Similarly, for a 3D CESE scheme using an unstructured spatial mesh built from tetrahedra, it can 
be shown that generally the time marching cannot be carried out independently over two disjoint sets of 
space-time staggered solution points. 

According to the above discussions, a CESE multidimensional time marching using an unstructured 
triangular or tetrahedral spatial mesh generally will involve all space-time solution points. This coupled 
with the fact that a CESE scheme has a space-time staggered stencil implies that solution decoupling may 
become a real nuisance. To deal with this problem once and for all, recently new CESE schemes have been 
developed such that possibility of solution decoupling can be removed completely. The description of these 
new schemes is one of the topics covered in the current paper. As will be shown, (i) these new schemes are 
constructed without compromising the basic CESE ideas, and (ii) their effectiveness has been verified by 
real-world applications. [62, 64-67]. 

To describe another topic of the current paper, note that the spatial solution point G ' depicted in Fig. 3 
is the centroid of the hexagon ABCDEF where points A , C and E, respectively, are rather arbitrarily 
designated as the centroids of the three neighbors of A FBD. In fact, to insure that CEs will not overlap in 
space and thus no loss of local and global flux conservation will occur, it only requires that points such as 
G, A, C and E lie in the interiors of A FBD, A FHB, ABLE, and A FDM, respectively. 

For the rather benign triangular domain depicted in Fig. 3, the hexagon ABCDEF happens to be 
convex, i.e., a line segment joining any two points in this hexagon always lies entirely in it. Because the 
centroid of a convex polygon always lies within its interior, the centroid (i.e., point G 1 ) of the hexagon 
ABCDEF lies within its interior. However, this may not be true for a more pathological triangular domain 
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formed from triangles with very large aspect ratios and sharply different orientations. As an example, 
consider Fig. 4 where A FBD is again surrounded by three neighboring triangles A FHB, A BLD and 
A FDM. However, in Fig. 4, the aspect ratio of each of the four triangles depicted is much larger than 
that of its counterpart depicted in Fig. 3. Moreover, in Fig. 4, the orientations of the sides BD and FD of 
A FBD are sharply different from those of the sides BH and FH of A FHB, respectively. Let points G, A, 
C and E depicted in Fig. 4 again be the centroids of A FBD, A FHB A BLD, and A FDM, respectively. 
Then the resulting hexagon ABCDEF may become highly concave and thus its centroid G' (not shown in 
Fig. 4) may no longer lie in its interior. 

Similarly, a CESE tetrahedral mesh is built from tetrahedra where any two neighboring tetrahedra 
share a common triangular face. Given any two neighboring tetrahedra, in the original CESE tetrahedral- 
mesh construction, the 5-vertex polyhedron with its vertices being (i) the centroids of these two neighbors 
and (ii) the three vertices of the triangular face shared by them, becomes the spatial projection of a BCE. 
Because each tetrahedron is associated with four neighbors, it is associated with four 5- vertex polyhedra and 
their union — a 8-vertex polyhedron. The 8-vertex polyhedron and its centroid, respectively, are the spatial 
projections of a CCE and a spatial solution point. For some pathological case where the 8- vertex polyhedron 
becomes highly concave, its centroid (i.e., the spatial solution point associated with this polyhedron) may 
lie outside of its interior. 

In a real-world CFD simulation near a solid wall, one often encounters a case where a boundary with 
high curvature or sharp corner is surrounded by triangular/tetrahedral meshes of extremely high aspect 
ratio (up to 10 6 ). For such an extreme case, a spatial solution point may lie far outside of the associated 
hexagon/8- vertex polyhedron. It could even be embedded beyond a solid-wall boundary and causes serious 
numerical difficulties. To overcome this problem, a new way of building spatial meshes from triangles and 
tetrahedra has been developed. Specifically, the centroids of triangles depicted in Fig. 3 have been replaced 
by the corresponding in-centers or other specially-defined interior points in the new construction. However, 
the solution point is still the centroid of the new hexagon formed by vertices and in-centers (or other 
specially-defined interior points). Similarly, centroids of tetrahedra also are replaced by corresponding in- 
centers or other specially-defined interior points in the new construction. According to the numerical results 
obtained [62,64,65], the new construction has the following advantages: (i) it can properly model curved 
boundaries — without it, computed surface gradient such as heat flux will be much less accurate (and noisier 
in many cases); (ii) it can reduce aspect ratios of triangular or tetrahedral cells used in gradient calculations 
in a highly stretched mesh and, as a result, can improve the overall accuracy for triangular and tetrahedral 
meshes; (iii) for vortex-shedding problems, better resolution is observed by using the in-center formulation. 
In this paper we will describe this new development and also provide the theoretical justification for the new 
construction. 

As mentioned earlier, mesh values of the physical dependent variables and their spatial derivatives are 
treated as independent mesh variables in a CESE scheme. As it turns out, for a case where a triangular or 
tetrahedral mesh is used, numerical evaluation of the spatial derivatives using the established CESE approach 
tends to become less accurate as the mesh cell aspect ratio becomes larger. In the latest development, the 
root cause behind this tendency was clearly identified and several remedies to overcome it were found through 
a rigorous mathematical analysis. However, because of the length of the current paper and the complexity 
of the mathematics involved, this latest work will be presented in another paper. 

Because it is rather simple and straightforward, in the following, first we will describe the new way of 
building spatial meshes from triangles and tetrahedra. The new flux-conserving solution-coupling procedures 
will be presented in Sec. 3. 

2. New construction of triangular and tetrahedral meshes 
2.1. Triangular meshes 

Again consider the triangle A FBD and the hexagon ABCDEF depicted in Fig. 4. As noted earlier, 
(i) because it is the centroid of the hexagon ABCDEF, the spatial solution point G' (not shown in Fig. 4) 
always lies in the interior of the hexagon if it is convex but may lie outside of it if it is concave; and (ii) 
serious numerical difficulties may arise if point G' lies outside of the hexagon ABCDEF. As such, to avoid 
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numerical difficulties, it is essential to ensure that ABC DEF is a convex hexagon. 

Because any triangle such as A FBD is always convex, the hexagon ABC DEF will also be convex if its 
boundary more or less coincides with that of A FBD. In other words, the hexagon is convex if points A, C, 
and E are located close enough to the sides BF, BD and FD, respectively. Unfortunately, because (i) A is 
the centroid of AFHB, and (ii) 


\BF\ < \BH\ and \BF\ < \FH\ (2.1) 

(hereafter, as an example, the length of BF is denoted by \BF |), point A generally is far away from the 
side BF. As such the hexagon ABC DEF may become concave. For this reason and the previously stated 
requirement that point A, C, and E must lie in the interiors of A FHB, A BLD, and A FDM, respectively, 
in the following, we will describe a new approach in which the locations of A, C, and E are chosen so that: 

(i) each of them again is an interior point of the associated triangle, and (ii) a point such as A will remain 
close to the side BF even for the pathological case Eq. (2.1). 

Consider Fig. 5. Here (i) points Pi, P 2 , P 3 are the vertices of the triangle APiP 2 P 3 which lies on the 
x-y plane and has its area 

A(AP 1 P 2 P 3 ) > 0 ( 2 . 2 ) 

(ii) Po is a point on the x-y plane to be defined shortly; (iii) for each k = 0, 1, 2, 3, Xk and yk are the x- and 
y- coordinates of point P&, respectively; (iv) point Q 1 is the point on P 2 P 3 (or its extension) which meets 
the condition 

TWi A TWs (2.3) 

i.e., Q 1 is the perpendicular projection of P 0 on the straight line which contains P 2 P 3 ; (v) point Q 2 is the 

point on P 3 P 1 (or its extension) which meets the condition 

PoQ 2 A JWi (2.4) 

i.e., Q 2 is the perpendicular projection of Po on the straight line which contains P 3 P 1 ; (vi) point Q 3 is the 

point on P 1 P 2 (or its extension) which meets the condition 

TWz A T\P ~2 (2.5) 

i.e., Q 3 is the perpendicular projection of Po on the straight line which contains P 1 P 2 ; (vii) a\, 0 . 2 , and <33 

are the internal angles of APiP 2 P 3 facing the sides P 2 P 3 , P 3 P 1 , and P 1 P 2 , respectively; (viii) 


def 


1 = |p 2 p 3 |, 


? 2 d =|P 3 Pi| and £3 d = |PiP 2 | 


(ix) 


def 


h k = \PoQk\ k = 1,2,3 

and (x) the x-y-z coordinate system depicted in Fig. 5 is a right-handed system, i.e., 

g x x Cy — e z , Gy x e z — g x and e z x g x Gy 


( 2 . 6 ) 

(2.7) 


(2.8) 


where e x , e y , and e z are the unit vectors in the x-, y-, and 2 - directions, respectively. 

Note that, according to Eq. (2.6) and Fig. 5, for each k = 1,2,3, £k denotes the length of the side of 
AP 1 P 2 P 3 facing the vertex P^. Also from the definition of points Q 1 , Q 2 , and Q 3 , and that of hk, k = 1,2, 3, 
one concludes that hi, / 12 , and h 3 , respectively, are the (shortest) distances from P 0 to the straight lines that 
contain the sides P2P3, P3P1, and P1P2, respectively. 

For each k = 0, 1, 2, 3, let ~P k be the position vector of point Pk, i.e., 




def _* _» 

— %k Vk 


k = 0, 1,2,3 


(2.9) 
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(2.10) 


Also let P o be an weighted average of P k, k = 1,2,3, i.e., 

~P 0 = W\P i + w 2 P 2 + w 3 P 3 

where w 1, W 2 , and w 3 are real weight factors with 

Wi+W2+W 3 = l (2-11) 

Moreover, given a pair of integers k-\ and fc 2 with k\ ^ fc 2 and fci,fc 2 = 0, 1,2,3, let Pk ± Pk 2 be the 
displacement vector joining points P/ C| and Pfc 2 and pointing in the direction from P^ to Pfc 2 , i.e., 

PkJ 5 k 2 = i ^k 2 -^k 1 , ki^k 2 and ki, k 2 = 0, 1, 2, 3 (2.12) 

Then, by using Eq. (2.12) and other well-known properties of the vector product such as 

~A x ~A = 0 and ~A x P = — P x Al (2.13) 

for any vectors A and P , one can show that 

P^ x P^ = (1*2 - l*i) X (i* 3 - Pi) = P 1 xP 2 + P 2 xP 3 + P 3 xP 1 (2.14) 
Similarly, with the aid of Eqs. (2.10)-(2.14), one can show that 

PPt X Pp ? 3 = Wi (Pi xP 2 + P 2 xP 3 + P 3 x 7*3) = Wi (pj ?2 X 1\P 3 ) (2.15) 

P0P3 X P0P1 = W2 (P 1 Xp 2 + P 2 Xp 3 + P 3 Xp 1^ = 1U 2 (P\P 2 x PiP 3 ) (2.16) 

P0P1 X P0P2 = »3 (?iX ?2 + ?2X?3 + ?3X = W 3 (p(P^ X P[P 3 ) (2.17) 

Next, by using (i) a well-known property of the vector product, (ii) the fact that A(APiP 2 Ps) denotes 

the area of APiP 2 P 3 , and (iii) Fig. 5, one has 

A(APiP 2 P 3 ) = (|^||^|sina 1 )/2= |P^ x P[P 3 1/2 (2.18) 

Similarly, with the aid of Eqs. (2.15)-(2.18), one can show that 

A(AP 0 P 2 P 3 ) = I Pp ?2 X PoH 1/2 = K| x A{APiP 2 P 3 ) (2.19) 

A(AP 0 P 3 Pi) = |iv£ x i^P[|/2 = |u> 2 | x A(AP!P 2 P 3 ) (2.20) 

and 

A{AP 0 PiP 2 ) = \PoP x ivt|/2 = |u, 3 | x A(AP!P 2 P 3 ) (2.21) 

where A(AP 0 P 2 P3), A(AP 0 P 3 Pi ), and A(AP 0 PiP 2 ) denote the areas of AP 0 P 2 P3, AP 0 P 3 Pi , and AP 0 PiP 2 , 
respectively. In turn, Eqs. (2.19)-(2.21) => 

A(AP 0 P 2 P 3 ) + A(AP 0 P 3 Pi) + A(AP 0 PiP 2 ) = (Kl + |w 2 | + I«7 3 |) X A(APiP 2 P 3 ) (2.22) 

Because a triangle such as A PiP 2 P 3 is intrinsically convex, according to a theorem on convex sets, point 
Po defined by Eqs. (2.10) and (2.11) lies on the boundary or in the interior of APlP 2 P3 -o- 

0 < tUfc < 1 k = 1,2,3 (2.23) 
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Moreover, because Eq. (2.2) implies that the interior of AP 3 P 2 P 3 is not empty, Po lies in the interior of 
AP 1 P 2 -P 3 

0 < < 1 A: =1,2,3 (2.24) 

For the case where Po lies on the boundary or in the interior of AP 1 P 2 P 3 , Eqs. (2.11), (2.22), and (2.23) 
now => the following obvious result: 

A(AP 0 P 2 P 3 ) + A(AP 0 P 3 Pi) + A(AP 0 PiP 2 ) = ^(AP!P 2 P 3 ) (2.25) 

On the other hand, for the case where Po lies outside of AP 1 P 2 P 3 , Eq. (2.23) is no longer true, i.e., some 
of the weight factors must be negative. Because Eq. (2.11) requires that some of the weight factors must be 
positive, we arrive at the conclusion that 

|wi| + \w 2 \ + \w 3 \ > \wi + w 2 + w 3 | = 1 (2.26) 

if Pq lies outside of APiP 2 P 3 . In turn, Eqs. (2.22) and (2.26) => 

A(AP 0 P 2 P 3 ) + A(AP 0 P 3 P 1 ) + A(AP 0 P 1 P 2 ) > AiAP^Ps) (2.27) 

if Pq lies outside of APiP 2 P 3 . 

Note that, for any real parameter (, Eqs. (2.2) and (2.6) => 

i k > 0 and (f k ) C ’ >0 k = 1,2,3; —00 < C < +00 (2.28) 

As such the weight factors 


= (yyrww * = 1 - 2 ' 3; -“<<<+“ < 2 - 29 > 

are well defined and satisfy the conditions 


0 < w k ( C) <1 k = 1,2,3, 

—00 < c < +°° 

(2.30) 

and 



wi( C) + w^ 2 (C) + w 3 (C) = 1 

— 00 < C < +00 

(2.31) 

In this subsection hereafter we consider only the special case where 


w k = w k ( C) A = 1 , 2 , 3; 

—00 < c < +°° 

(2.32) 


For this special case, Eqs. (2.30)-(2.32) imply that both Eqs. (2.11) and (2.24) are satisfied for any real 
number £. Thus, for the case Eq. (2.32), (i) point Pq defined by Eq. (2.10) must always lie in the interior of 
AP 1 P 2 P 3 , and (ii) Eqs. (2.19)-(2.22) reduce to 


A(APoP 2 P 3 ) — w 1 x A(APiP 2 P 3 ) 


(2.33) 

A(APoP 3 P\) = w 2 x A(APiP 2 P 3 ) 


(2.34) 

A(APoPiP 2 ) = w 3 x A(AP\P 2 P 3 ) 

and Eq. (2.25), respectively. Moreover, for the special case £ = 0, we have 


(2.35) 

Wi = w 2 = w 3 = 1/3 and 1^ 0 = (~^ 1 + ^2 + ^ 3 ) /3 

(C = 0 ) 

(2.36) 
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Because the position vector of the centroid of a triangle is the simple average of those of its three vertices, 
Eq. (2.36) => point Po is the centroid of APlP 2 P 3 if £ = 0. 

Note that, for a polygon other than a triangle (which is the simplest polygon), generally it is not true 
that the position vector of its centroid is the simple average of those of its vertices. As an example, the 
position vector of the centroid of a quadrilateral may not be the simple average of those of its four vertices 
unless the quadrilateral is a parallelogram. 

Next, with the aid of Eqs. (2.6) and (2.7), and Fig. 5, one can see that 

A(AP 0 P 2 P 3 )=^1, A(AP 0 P 3 Pi)=^ and A(AP 0 P!P 2 ) = ^ (2.37) 

In turn, by substituting Eq. (2.37) into Eqs. (2.33)-(2.35) and using Eqs. (2.29) and (2.32), one has 


= w k A(AP 1 P 2 P 3 ) 


(4) c x A(AP 1 P 2 P 3 ) 

(ttf + (4)< + (4) c 


k = 1,2,3; —00 < £ < + 00 (2.38) 


and thus 

hi h- 2 h 3 2 x A(APiP 2 P 3 ) 

~ ~ (^F 1 “ (4F + (4)< + (4) c 


— 00 < £ < +00 (2.39) 


Because (i) h 1 , h 2 and h 3 are the distances from point Po to the sides P 2 P 3 , P 3 Pi, and PiP 2 , respectively, 
and (ii) 4, 4, and £3 are the lengths of P 2 P 3 , P 3 Pi, and P 3 P 2 , respectively, Eq. (2.39) states that, for a 
given APiP 2 P 3 and a given £, the distance from point Pq to any side of APiP 2 P 3 is proportional to the 
(£ — 1 )th power of the length of this side. 

In fact, for the case £ = 0 (i.e., when P 0 is the centroid of APiP 2 P 3 ), Eq. (2.39) reduces to 


hk 


'2 x A{APiP 2 P 3 )' 


1 

£k 


k = 1,2,3 if £ = 0 


(2.40) 


i.e., the distance between the centroid Pq and a side of A P 3 P 2 P 3 is inversely proportional to the side length. 
Thus the centroid is farthest away from the shortest side and closest to the longest side — a result consistent 
with what we observe from the triangle A FHB depicted in Fig. 4. 

On the other hand, for the case £ = 1, Eq. (2.39) reduces to 


hi — h 2 — h 3 


2 x A{APiP 2 P 3 ) 
£\ + £2 + £ 3 


if £ = 1 


(2.41) 


i.e., point Po has the same distance to the three sides of A P 3 P 2 P 3 if £ = 1. By definition, this implies that 
Po is the in-center of A PlP 2 P 3 for the case £ = 1. 

In general one can infer from Eq. (2.39) that: 

(a) For the case £ > 1, point Po is farthest away from the longest side of A P 3 P 2 P 3 and closest to its shortest 
side. The trend becomes stronger as the value of £ increases from 1. 

(b) For the case £ < 1, point P 0 is farthest away from the shortest side of A PiP 2 P 3 and closest to its longest 
side. The trend becomes stronger as the value of £ decreases from 1. 

Let the interior point Po of any triangle APiP 2 P 3 defined by Eqs. (2.10) and (2.32) be referred to as 
the interior point of this triangle associated with the parameter £. Then, from the above observations and 
the discussions given in Sec. 1, it becomes clear that a spatial solution point such as the point G' depicted 
in Fig. 3 will be much less likely to lie outside of the associated hexagon ABC DEF and causes numerical 
difficulties if each of the centroidal mesh points A, C, G and E depicted in Figs. 3 and 4 is replaced by a an 
interior point associated with a parameter £ > 1. In fact, it has been shown through numerical experiments 
that the numerical difficulties associated with spatial solution points being embedded beyond a solid-wall 
boundary can be overcome by replacing centroids of triangles with in-centers (i.e., interior points associated 
with £ = 1) in the construction of a spa tial triangular mesh. 
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2.2. Tetrahedral meshes 

Consider a 3D spatial domain formed by tetrahedra with the understanding that any two neighboring 
tetrahedra share a common triangular face. As an example, consider the tetrahedron P1P2P3P4 depicted in 
Fig. 6(a) where Pi, P2, P3 and P4 are the vertices of the tetrahedron. The tetrahedron has four triangular 
faces, i.e., AP1P2P3, AP2P3P4, AP3P4P1, and AP4P1P2. As shown in Fig. 6(b), the tetrahedron P1P2P3P4 
shares these faces with four neighboring tetrahedra, i.e., P1P2P3D, P2P3P4A, P3P4P1B, and P4P1P2C, 
respectively. Note that, to simplify notation, the same symbols used in Sec. 2.1 may also be used in this 
subsection albeit that those used here generally denote objects in the 3D x-y-z space. 

Because a tetrahedron is intrinsically convex in a 3D x-y-z space, its centroid must lie in its interior if 
its volume > 0. Hereafter, only tetrahedra of volume > 0 will be considered. In particular, we assume that 

Hpf 

V (P\ P2P3P4) = the volume of the tetrahedron P1P2P3P4 > 0 (2.42) 

Let the points Po, A 0 , B 0 , Co, and D 0 depicted in Fig. 6(b) be the centroids of the tetrahedra, P1P2P3P4, 
P2P3P4A, P3P4P1P, P4P1P2C, and P1P2P3P, respectively. Then each of these five centroids lies in the 
interior of the associated tetrahedron. As such the four polyhedra P0P2P3P4A0, P0P3P4P1B0, P0P4P1P2C0, 
and P0P1P2P3P0 (each of them has five vertices, six triangular faces, and nine edges) do not overlap. In 
the original CESE 3D tetrahedral mesh construction, (i) each of these 5-vertex polyhedra is the spatial 
projection of a BCE, (ii) the union of the above four 5-vertex polyhedra is a polyhedron with 8 vertices, 
12 triangular faces, and 18 edges, and also is the spatial projection of a CCE, and (iii) the centroid of the 
8-vertex polyhedron referred to above is a spatial solution point, i.e., the spatial projection of a space-time 
solution point. 

As in the 2D case, for some pathological case where an 8-vertex polyhedron become highly concave, its 
centroid (i.e., the associated spatial solution point) may lie outside of the polyhedron and causes numerical 
difficulties. In the following, we explain how these difficulties can be overcome by a new construction of 3D 
tetrahedral meshes. 

Consider Fig. 6(a) again. Here (i) Po is a point in the x-y-z space to be defined shortly and not necessary 
the centroid of the tetrahedron P1P2P3P4; (ii) for each k = 0, 1, 2, 3,4, Xk, y and z i~ are the x-, y- and ,2- 
coordinates of point Pfc, respectively; (iii) point Q± (not shown in Fig. 6(a)) is the point that lies on the 
plane containing the face AP2P3P4, and also meets the condition 


P0Q1 1 AP 2 P 3 P 4 (2.43) 

i.e., Q 1 is the perpendicular projection of Po on the plane containing AP2P3P4; (iv) point Q 2 (not shown) 
is the point that lies on the plane containing the face AP3P4P1, and also meets the condition 


P0O2 JL AP3P4P1 (2.44) 

i.e., Q 2 is the perpendicular projection of Po on the plane containing AP3P4P1; (v) point Q3 (not shown) is 
the point that lies on the plane containing the face AP4P1P2, and also meets the condition 


P0Q3 1 AP 4 P 1 P 2 (2.45) 

i.e., Q 3 is the perpendicular projection of Po on the plane containing AP4P1P2; (vi) point Q4 (not shown) 
is the point that lies on the plane containing the face AP1P2P3, and also meets the condition 

Wh A AP4P 2 P3 (2.46) 

i.e., Q4 is the perpendicular projection of Po on the plane containing AP1P2P3; (vii) 

A 1 = A(AP 2 P 3 P 4 ), A 2 d = A(AP 3 P 4 P 1 ), A 3 d A f A(AP 4 PiP 2 ), and A 4 = AiAP^Ps) (2.47) 
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where A(AP 2 P 3 P 4 ), A(AP 3 P4 Pi), A(AP 4 PiP 2 ), and A(APiP 2 P 3 ) are the areas of the four triangular faces 
AP2P3P4, AP3P4P1, AP4P1P2, and AP1P2P3, respectively; (viii) 


h k = \PoQk\ k = 1,2, 3, 4 (2.48) 

and (ix) the x-y-z coordinate system depicted in Fig. 6(a) is also a right-handed system which satisfies 
Eq. (2.8). 

Note that, according to Eq. (2.47) and Fig. 6(a), for each k = 1,2, 3, 4, A k denotes the area of the 
face of the tetrahedron P1P2P3P4 which is opposite to the vertex P k . Also from the definition of points 
Q 1, Q 2, Q 3, and Q 4 and that of h k , k = 1,2, 3, 4, one concludes that h\, h 2 , J13, and /14, respectively, are 
the (shortest) distances from Pq to the planes that contain the faces AP2P3P4, AP3P4P1, AP4P1P2, and 
AP1P2P3, respectively. 

For each k = 0, 1,2, 3, 4, let P k be the position vector of point Pk, i.e., 

P k = f x k e x + y k e y + z k e z A: = 0,1, 2, 3, 4 (2.49) 

Also let P 0 be an weighted average of P k , k = 1,2, 3, 4, i.e., 

P 0 = w{P \ + w 2 P 2 + w 3 P 3 + W 4 P 4 (2.50) 

where w\, w 2 , w 3 , and W 4 are real weight factors with 

w\ + w 2 + w 3 + W 4 = 1 (2-51) 

Moreover, given a pair of integers k\ and k 2 with k\ ^ k 2 and k -\ , k 2 = 0, 1,2, 3, 4, let Pk 1 P k2 be the 
displacement vector joining points P kl and P k2 and pointing in the direction from P k 1 to Pfc 2 , i.e., 

P^K 2 d =V k2 -^ kl , h^k 2 and k u k 2 = 0,1, 2, 3, 4 (2.52) 


Then, by using Eq. (2.52) and other well-known properties of the vector product such as Eq. (2.13), one can 
show that Eq. (2.14) is still valid in the current 3D case. In turn, with the aid of Eqs. (2.14) and (2.52), and 
other well-known properties of the scalar triple product such as 

(1 x If) • ~6 = 0 x (?) • ~A = (Cl x ~1) ■ T$ = -(1 x7J)-~8 (2.53) 

and 

(1 x 'l) ■ H = (1 x ~§) ■ = 0 x It) • = 0 (2.54) 

for any vectors , and , one can show that 


A = f (Pj? 2 X P^) = 4 X ~P 2 ) ■ ~P 4 + (P 2 X P3) ■?4 + (?3X? 1 )'?4 + (P 3 X P 2) • P 1 (2.55) 


Moreover, by using Eqs. (2.13), (2.50) and (2.52), one has 

P 0 Pl X P 0 P 2 = (1 - Wi - W 2 )P i xP 2 + ( W3P3 + wiPij x (Pi - P 
Next, by using Eqs. (2.50)-(2.56), one can show that 


(2.56) 


(P 0 Pi x P 0 P 2 ) • P 0 P 3 


= W4 


(P4 X P 2 ) • P 3 + (P 2 X P4) • P 4 + (P 3 X P 4) • P 1 + (1* 4 X 1*3) ' P-2 


= — U>4 A 


(2.57) 


N ASA/TM— 20 1 3-218061 


15 



In turn, by carrying out the index mapping 


(1j 2,3,4) — > (2,3,4, 1) 
of Eq. (2.57) repeatedly and using Eqs. (2.53)-(2.55), one obtains 




{P 0 P 2 x P 0 P 3 ) • P 0 P 4 


= W 1 


(?2 X 1*3) ■ ^4 + (1*3 X P 2 ) ■ i*l + X P 4 ) • P 2 + (l*, X P 4 ) ■ Po 


» ~ , ~ ~ > 


(P 0 P 3 X P 0 P 4 ) • P 0 Pl 


= w 2 


(P 3 X P 4 ) ■ Pi + ( 1 * 4 x P 3) ■ P 2 + {P 1 xP 2 )-P 3 + (P 2 X ~P 1 ) • ~P 4 


= w\ A 


= —w 2 A 


and 


(Jvt x ivt) • pH 
= ^3 [(^ 4 X Pi) • P 2 + (P 1 xP 4 )-P 3 + (P 2 x 7 * 3 ) • P4 + (P 3 X P 2 ) * P 1 


= W 3 A 


(2.58) 


(2.59) 


(2.60) 


For any four different integers k 4 ,k 2 , k 3 , k 4 = 0, 1, 2, 3, 4, let 

def 

V ( Pfc, P/ C2 Pfe 3 Pk 4 ) = the volume of the tetrahedron with the vertices Pfcj , Pfc 2 , Pfc 3 , and P/ C4 (2.61) 
Then by using the definition of the scalar triple product along with Eqs. (2.55) and (2.57)-(2.60), one has 

V(P\P 2 P 3 P 4 ) = | (1\P 2 x I\P 3 ) ■ P[P 4 1/6 = |A|/6 (2.62) 

V(P 0 P 2 P 3 P 4 ) = \(1PP 2 x Ppp) ■ P^P 4 1/6 = \w 4 A|/6 = \w 4 \ x V (P\P 2 P 3 P 4 ) (2.63) 

V{PqP 3 P 4 P\) = |(iv^ x ivt) • FVPlI/6 = \w 2 A |/6 = |w 2 | x V{P 4 P 2 P 3 P 4 ) (2.64) 

V(P 0 P 4 PiP 2 ) = \(PoP 4 x ivt) • iV^|/6 = \w 3 A |/6 = \w 3 \ x I/(PiP 2 P 3 P 4 ) (2.65) 

and 

V(P 0 PiP 2 P 3 ) = I (IpPi X ivt) • PoP\/6 = \w 4 A | / 6 = |zc 4 | x V(P 4 P 2 P 3 P 4 ) (2.66) 

Eqs. (2.63)-(2.66) =4- 

V(P 0 P 2 P 3 P 4 ) + E(P 0 P 3 P 4 Pi) + V(P 0 P 4 PiP 2 ) + V(P 0 PiP 2 P 3 ) 

= (|wh| + |?U 2 | + 1^3 1 + 1^4 1 ) X y^PlPpPi) 

Because the tetrahedron P 4 P 2 P 3 P 4 is intrinsically convex in the x-y-z space, point Pq defined by 
Eqs. (2.50) and (2.51) lies on the boundary or in the interior of the tetrahedron P 1 P 2 P 3 P 4 4=> 

0 < w k < 1 k = 1,2, 3, 4 (2.68) 

Moreover, because Eq. (2.42) implies that the interior of the tetrahedron P 4 P 2 P 3 P 4 is not empty, Pq lies in 
the interior of the tetrahedron P\P 2 P 3 P 4 <t=> 

0<w k <l k = 1,2, 3, 4 (2.69) 

For the case where Po lies on the boundary or in the interior of the tetrahedron P 4 P 2 P 3 P 4 , Eqs. (2.51), 
(2.67), and (2.68) now => the following obvious result: 

V(P 0 P 2 P 3 P 4 ) + y(P 0 P 3 P 4 Pi) + E(P 0 P 4 PiP 2 ) + V(P 0 PiP 2 P 3 ) = V(P 4 P 2 P 3 P 4 ) (2.70) 
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On the other hand, for the case where Pq lies outside of the tetrahedron P 1 P 2 P 3 P 4 , Eq. (2.68) is no longer 
true, i.e., some of the weight factors must be negative. Because Eq. (2.51) requires that some of the weight 
factors must be positive, we arrive at the conclusion that 

|wi| + \w 2 | + 1 103 1 + \Wi\ > |u?i + W 2 + W 3 + 104 1 = 1 (2-71) 

if Pq lies outside of the tetrahedron P 1 P 2 P 3 P 4 . In turn, Eqs. (2.67) and (2.71) => 


V(P 0 P 2 P 3 P 4 ) + y(PoP 3 P 4 Pi) + E(P 0 P 4 PiP 2 ) + E(P 0 PiP 2 P 3 ) > E(PiP 2 P 3 P 4 ) (2.72) 


if Pq lies outside of the tetrahedron P 1 P 2 P 3 P 4 . 

Note that, for any real parameter £, Eqs. (2.42) and (2.47) 


Ak > 0 and (A*,)^ > 0 
As such the weight factors 

(4fc) C 


k = 1 , 2 , 3, 4; —00 < C < +00 


(2.73) 


t def 

Wk{Q = 


(Ai)^ + (A 2 )C + (A 3 )C + (^. 4 )? 

are well defined and satisfy the conditions 

0 < Wk(C) < 1 fc = l,2,3,4 — 00 < ( < +00 


fc = l,2,3,4; — 00 < ( < +00 (2.74) 


(2.75) 


and 

«d(0 + w 2 (C) + tU3(C) + ^4(C) = 1 -oo<C<+oo (2.76) 

Note that, even though the same notation is used, the weight factors Wk(C )> k = 1, 2, 3, 4, defined in Eq. (2. 74) 
are different from those defined in Eq. (2.29). The definition Eq. (2.74) will be used in this subsection. 
Moreover, hereafter we consider only the special case where 

Wk = Wk(0 k= 1,2, 3,4; — 00 < £ < +00 (2.77) 

For this special case, Eqs. (2.75)-(2.77) imply that both Eqs. (2.51) and (2.69) are satisfied for any real 

number £. Thus, for the case Eq. (2.77), (i) point Pq defined by Eq. (2.50) must always lie in the interior of 

the tetrahedron P 1 P 2 P 3 P 4 , and (ii) Eqs. (2.63)-(2.67) reduce to 


V(P 0 P 2 P 3 P 4 ) =w 1 x V ( P| P 2 P 3 P 4 ) 


(2.78) 

V(PqP 3 P a Pi) =w 2 x V ( P| P 2 P 3 P 4 ) 


(2.79) 

E(P 0 P 4 PiP 2 ) = w 3 x E(PiP 2 P 3 P 4 ) 


(2.80) 

V (PqPiP 2 P 3 ) =w 4 x E(P 1 P 2 P 3 P 4 ) 

and Eq. (2.70), respectively. Moreover, for the special case £ = 0, we have 


(2.81) 

uq = w 2 = w 3 = w 4 = 1/4 and ~P 0 = + ~P 2 + ~P 3 + /4 

(C = o) 

(2.82) 


Because the position vector of the centroid of a tetrahedron is the simple average of those of its four vertices, 
Eq. (2.82) =>■ point Pq is the centroid of the tetrahedron P 1 P 2 P 3 P 4 if £ = 0. 

Note that, for a polyhedron other than a tetrahedron (which is the simplest polyhedron), generally it is 
not true that the position vector of its centroid is the simple average of those of its vertices. As an example, 
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the position vector of the centroid of a hexahedron may not be the simple average of those of its eight vertices 
unless the hexahedron is a parallelepiped. 

Next, with the aid of Eq. (2.47) and (2.48), and Fig. 6 (a), one can see that 

E(P 0 P 2 P 3 P 4 ) = V(P 0 P 3 P 4 Pi) = E(P 0 P 4 PiP 2 ) = ^ and V(P 0 PiP 2 P 3 ) = (2.83) 

In turn, by substituting Eq. (2.83) into Eqs. (2.78)-(2.81) and using Eqs. (2.74) and (2.77), one has 


h f k = w fe y(PiP 2 p 3 P4) 


(A k )i x E(P 1 P 2 P 3 P 4 ) 
(-^i)^ + (A2Y + (.A 3 )£ + (A 4 ) < ’ 


fc = l,2,3,4; — 00 < £ < +00 (2.84) 


and thus 

hi _ h 2 _ h 3 _ hi _ 3 x V(PiP 2 P 3 P 4 ) 

_ (A^T - (ZOC^T - (A,)C + (A 2 )C + (A 3 )C + (A 4 )< 


— 00 < £ < +00 (2.85) 


Because (i) hi, h 2 , h 3 , and hi are the distances from point Po to the triangular faces A P 2 P 3 P 4 , A P 3 P 4 P 4 , 
A P 4 P 1 P 2 , and A P\P 2 P 3 respectively, and (ii) Ai, A 2 , A 3l and A 4 are the areas of A P 2 P 3 P 4 , A P 3 P 4 P 1 , 
AP 4 PiP 2 , and APiP 2 P 3 , respectively, Eq. (2.85) states that, for a given tetrahedron P\P 2 P 3 P 4 and a given 
£, the distance from point Pq to any face of this tetrahedron is proportional to the (£ — 1 )th power of the 
area of the face. 

In fact, for the case £ = 0 (i.e., when P 0 is the centroid of the tetrahedron PiP 2 P 3 P 4 ), Eq. (2.85) reduces 
to 

, _ [3 x E(PiP 2 P 3 P 4 ) 


1 

A k 


k= 1,2, 3, 4, if £ = 0 


(2.86) 


i.e., the distance between the centroid Pq and a face of the tetrahedron PlP 2 P 3 P 4 is inversely proportional 
to the area of this face. Thus the centroid is farthest away from the smallest face and closest to the largest 
face. 

On the other hand, for the case £ = 1, Eq. (2.85) reduces to 


hi — h 2 — h 3 — h 4 


3 x E(PiP 2 P 3 P 4 ) 
Ai + A 2 + ^3 + A 4 


if C = 1 


(2.87) 


i.e., point Po has the same distance to the four faces of the tetrahedron PiP 2 P 3 P 4 if ( = 1. By definition, 
this implies that Pq is the in-center of the tetrahedron PiP 2 P 3 P 4 for the case £ = 1 . 

In general one can infer from Eq. (2.85) that: 

(a) For the case £ > 1, point Po is farthest away from the largest face of the tetrahedron PiP 2 P 3 P 4 and 
closest to its smallest face. The trend becomes stronger as the value of ( increases from 1. 

(b) For the case ( < 1, point Pq is farthest away from the smallest face of the tetrahedron PiP 2 P 3 P 4 and 
closest to its largest face. The trend becomes stronger as the value of £ decreases from 1. 

Let the interior point Po of any tetrahedron PiP 2 P 3 P 4 defined by Eqs. (2.50) and (2.77) be referred to 
as the interior point of this tetrahedron associated with the parameter £. Then, from the above observations 
and the discussions given in Sec. 1, it becomes clear that a spatial solution point which is the centroid of 
the 8-vertex polyhedron referred to earlier will be much less likely to lie outside of the associated 8-vertex 
polyhedron and causes numerical difficulties if each of the centroidal mesh points Pq, Aq, Bq, Cq and Dq 
depicted in Fig. 6(b) is replaced by an interior point associated with a parameter £ > 1. In fact, it has been 
shown through numerical experiments that the numerical difficulties associated with spatial solution points 
being embedded beyond a solid-wall boundary can be overcome by replacing centroids of tetrahedra with 
in-centers (i.e., interior points associated with £ = 1 ) in the construction of a spatial tetrahedral mesh. 
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3. Flux-conserving solution-coupling procedures 

In this section, it will be described how an ID CESE scheme can be modified so that the two decou- 
pled solutions of the original scheme can become coupled in the new scheme while still preserving a flux 
conservation condition to be defined shortly. This description is then followed by a sketch of the 2D and 3D 
extensions. 

3.1. ID CESE solution-coupling schemes 

As a preliminary, note that the CEs, SEs and the ID schemes, which were defined over f2i in Sec. 1, 
hereafter will be defined over 


ft = {(j, n)\j, n = 0, ±1/2, ±1, ±3/2, . . .} = f2i U f2 2 (3.1) 

In particular, it should be understood that the symbol fli from now on is replaced by 12 in each of Eqs. (1.3)- 
(1.16). Also, as an example, the extension (defined over Q) of the a-e scheme will be referred to as the dual 
a-e scheme. 

To clarify the significance of the extension, we offer the following observations: 

(a) The solution of a dual scheme is formed by two decoupled solutions which are defined over 12 1 and D 2 , 
respectively. 

(b) The same space-time region can be designated as a BCE of a mesh point £ 12i and also as a BCE of 
another mesh point £ 122- For examples, the space-time region ABCD depicted in Fig. 2(c) can be 
designated as CE_(j,n) and CE + (j — 1/2, n), respectively, while the region AFED depicted in Fig. 2(d) 
as CE_| -(j,n) and CE_(j ± 1/2, n), respectively. Hereafter, as an example, points ( j,n ) and (j — 1/2, n) 
(i.e., points A and B) will be referred to as the hosts of the region ABCD. 

(c) Let the region ABCD be designated as CE_(j, n). Then, on its boundary, the line segments AB and 
AD belong to SE(j, n) while CB and CD belong to SE(j — 1/2, n — 1/2), i.e., over AB and AD, the 
vector h* should be evaluated using the solution data stored at point A while, over CB and CD, it 
should be evaluated using those stored at point C. 

Alternately, let the region ABCD be designated as CE + (j — 1/2, n). Then, on its boundary, the line 
segments BA and BC belong to SE(j — 1/2, n) while DA and DC belong to SE(j, n— 1/2), i.e., over BA 
and BC, the vector h* should be evaluated using the solution data stored at point B while, over DA 
and DC, it should be evaluated using those stored at point D. Note that, in the above and hereafter, 
as an example, a line segment joining points A and B is denoted by AB if it belongs to SE(A) (i.e., the 
SE centered at point A) while the same line segment is denoted by BA if it belongs SE(B) (i.e., the SE 
centered at point B ). 

(d) From the discussions given in the above item (c), it becomes clear that 

(f h* ■ ds = 0 and (f h* ■ ds = 0 (j, n) £ fi (3.2) 

/s(CE^(j,n)) is(CE_(j+l/ 2 ,n)) 

represent two different flux conservation conditions even though CE + (j, n) and CE_ (j ±1/2 ,n) occupy 
the same space-time region for any (j, n) £ 12. 

As will be shown, each marching step of the ID solution-coupling scheme is formed by two sub-steps. 
In the first sub-step, referred to as the dual-scheme marching step, the intermediate mesh values at the 
nth time level are evaluated from the given mesh values at the (n — l/2)th time level using a dual-scheme 
(see Fig. 7). In the second sub-step, referred to as the solution-coupling step, the final mesh values at the 
nth time level are updated from the intermediate mesh values through a flux-conserving weighted-averaging 
procedure. Let the intermediate mesh values of u and its normalized spatial derivative Ux at any (j, n) £ fi 
be denoted by h" and (h®)", respectively. Then, as an example, ft" and (-u s )" can be evaluated in terms of 
the given mesh values at the (n — l/2)th time level by using the dual a-e scheme defined by Eqs. (1-12) and 
(1.15), i.e., 

Uj = ^ { (! - v) [u - (! ± ^)u $]”+ 1 1 / / 2 ± (1 ± v) [« ± (1 - v)ux] n rl',t } ( j , n) £ 12 (3.3) 
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and 


(%)" = \ { K 1 - e)u + (2e - 1 - ^Us]" + i /2 “ K 1 - e)u + (1 - v - 2e)u®]”_ 1 1 / / 2 2 j (j, n ) G fi (3.4) 

Note that, by using a relation similar to the second equation in Eq. (1.9), one can define the “unnormalized” 
version of (it®)", he., 

(i.n)efi (3.5) 

J AX J 

To a void introducing unnecessary new notations, in the following description of the solution-coupling 
step, at each ( j , n) G O, the final mesh values after updating from the intermediate values are simply denoted 
as uj and (u s )", respectively. 

Note that the line segment AF depicted in Fig. 7 belongs to SE (j,ri). Thus, by using (i) Eqs. (1.4), 
(1.6) and (1.9), and (ii) Fig. 2(d), one concludes that the flux of h* leaving CE +(j,ri) through AF evaluated 
using the final solution data stored at point A is given by 

F(AF) = / u*(x,t n ;j,n)dx=—(u + )‘j (j,n)&Q (3.6) 

J Xj 

where 

(«+)" = f (u + u®)^ (j, n) G (3.7) 

Let F(AF) be the intermediate value of F(AF), i.e., the value of F(AF) when it" and (it®)" are replaced 
by v,™ and (u®)", respectively. Then Eqs. (3.5)-(3.7) => 

F(AF) = ^(u + )] (j.n)Gfi (3.8) 

where 

(w+)" = f (w + «®)" (i, n) G 0 (3.9) 

At this juncture, note that: (i) at the midpoint M + of AF, x = Xj + (ax/4) and t = t n \ and (ii) 
Eqs. (1.6), (1.9), and (3.7) => 

u*(xj + (ax/4 ),t n ;j,n) = (u+)" (j,n) G fi (3.10) 

Thus, (u+)" is the value of it at the midpoint M + of AF evaluated using the values it" and (it®)" with 
the aid of Eqs. (1.6) and (1.9). Similarly (u + )" is the value of u at the midpoint M + evaluated using the 
intermediate values u " and (it®)" with the aid of Eqs. (1.6), (1.9), and (3.5). 

Because (it + )" and (it+)" are evaluated using the solution data stored at point (j, n) (i.e., point A), 
hereafter, respectively, they will be referred to as the final and intermediate values of u at the midpoint 
M + associated with the host point A of the region AFED. Similarly, F(AF) and F(AF), respectively, are 
referred to as the final and intermediate values of the flux leaving the top face of the same region associated 
with the same host. 

Moreover, with the aid of Eqs. (3.6) and (3.8), and Fig. 7, one concludes that F(AF) ( F(AF )) is equal 
to the length of AF multiplied by the approximated value (u + )" ((ti + )") of u at the midpoint M + . In other 
words, the value of u at the midpoint M + is the average value of u over AF — a result stemmed from the 
fact that u*(x,t;j,n) is a first-order Taylor’s expansion in x when t = t n (see Eqs. (1.3) and (1.6)). 
Similarly, with the aid of Fig. 7, one concludes that: 

(a) The flux of h* leaving CE_ ( j, n ) through the line segment AB evaluated using the final solution data 
stored at point A is given by 

, t f x j — ( Ax/2) , 

F(AB) = / u*(x,t n \j,n) dx = —(«,_)" (j,n) G (3.11) 

J Xj 
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where 

(«-)? d = (u - u *)!? (j, n) £ O (3.12) 

is the value of u at the midpoint M_ of AB evaluated using the values it” and (it $ )” with the aid of 
Eqs. (1.6) and (1.9). 

Also F(AB), the intermediate value of F(AB), is given by 

F(AB) = (3.13) 

where 

(u_)7 d ^(«-fis)7 (j,n)eil (3.14) 

is the value ofu at the midpoint M_ evaluated using the intermediate values ft” and (h 2 )" with the aid 
of Eqs. (1.6), (1.9), and (3.5). 

Because (u_)” and (h_)" are evaluated using the solution data stored at point (j, n) (i.e., point A), 
hereafter, respectively, they will be referred to as the final and intermediate values of u at the midpoint 
M_ associated with the host point A of the region ABCD. Similarly, F(AB) and F(AB), respectively, 
are referred to as the final and intermediate values of the flux leaving the top face of the same region 
associated with the same host. 

Moreover, with the aid of Eqs. (3.11) and (3.13), one concludes that F(AB) (F(AB)) is equal to the 
length of AB multiplied by the approximated value (w_)” ((«_)”) of u at the midpoint M_. 

(b) The flux of h* leaving CE_(j + 1/2, n) through the line segment FA evaluated using the final solution 
data stored at point F is given by 

F(FA) = J u*(x,t n \j + 1/2, n)dx = — (w_)” +1/2 (j,n) £ O (3.15) 

where 

( u -)j+ 1 /2 d = ( u - «x)"+i/2 (i . n ) e n (3.16) 

is the value ofu at the midpoint M + of FA evaluated using the values w" +1 / 2 and (%)" +1 , 2 with the aid 
of Eqs. (1.6) and (1.9). Note that: (i) because (j,n) £ fl <t=> (j + 1/2, n) G O, Eq. (3.12) Eq. (3.16); 
and (ii) according to Eqs. (3.6) and (3.15), F(AF) ^ F(F A) even though AF and FA denote the same 
line segment. In fact, F(AF) is a flux evaluated using the solution data stored at point A while F(F A) 
is another flux evaluated using those stored at point F. The reader is warned that other similar pairs 
of notations will appear throughout the paper. 

Also F(F A), the intermediate value of F(F A), is given by 

^ At 

F(FA) = —(u.)] +1/2 (j,n)£ n (3.17) 

where 

(h-)" + 1 /2 = f (u - «*)"+i /2 (/. n) £ n (3.18) 

is the value ofu at the midpoint M+ evaluated using the intermediate values 'u” +1 , 2 and (w s )"+ 1/2 with 
the aid of Eqs. (1.6), (1.9), and (3.5). Obviously Eq. (3.14) <t=> Eq. (3.18). 

Because (w_)" +1 , 2 and (h_)” +1 ^ 9 are evaluated using the solution data stored at point (j + 1/2, n ) (i.e., 
point F), hereafter, respectively, they will be referred to as the final and intermediate values of u at the 
midpoint M + associated with the host point F of the region AFED. Similarly, F(FA) and F(F A), 
respectively, are referred to as the final and intermediate values of the flux leaving the top face of the 
same region associated with the same host. 
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Moreover, with the aid of Eqs. (3.15) and (3.17), one concludes that F(FA) ( F(F A )) is equal to the 
length of FA multiplied by the approximated value (u_)" +1 , 2 ((u_)” +1 , 2 ) of it at the midpoint M+. 

(c) The flux of h* leaving CE + (j — 1/2, n) through the line segment BA evaluated using the final solution 
data stored at point B is given by 

F{BA) = f f u*(x,t n \j -l/2,n)dx = ^(u+)]_ 1/2 (j,n)GO (3.19) 

Jxj-( Ax/2) 1 

where 

(«+)"_ 1/2 = f (u + «x)-_ 1/2 (j, n) G O (3.20) 

is the value of u at the midpoint M_ of RA evaluated using the values u r )j_ 1 / 2 and ( u x)”_i/ 2 with the 
aid of Eqs. (1.6) and (1.9). Obviously, Eq. (3.7) <t=> Eq. (3.20). 

Also F(BA), the intermediate value of F(BA), is given by 

^ At 

F(BA) = — (u+)?_ 1/2 (j, n) G n (3.21) 

where 

(u+)"_ 1/2 = f (n + u,)]_ 1/2 (j, n) G O (3.22) 

is the value of u at the midpoint M_ evaluated using the intermediate values u 1 l_ 1 / 2 and (%)"_ 1( / 2 with 
the aid of Eqs. (1.6), (1.9), and (3.5). Obviously Eq. (3.9) <t=> Eq. (3.22). 

Because (w + )"_ 1 y 2 and (h + )”_ 1 ^ 2 are evaluated using the solution data stored at point (j — 1/2, n) (i.e., 
point B), hereafter , respectively, they will be referred to as the final and intermediate values of u at the 
midpoint M_ associated with the host point B of the region ABCD. Similarly, F(BA) and F(BA), 
respectively, are referred to as the final and intermediate values of the flux leaving the top face of the 
same region associated with the same host. 

Moreover, with the aid of Eqs. (3.19) and (3.21), one concludes that F(BA) ( F(BA )) is equal to the 
length of BA multiplied by the approximated value (u+)/_ 1 / 2 ((tt+)"_ 1 / 2 ) of u at the midpoint M_. 
Next recall that two different final fluxes F(AF) and F(FA), respectively associated with the hosts 
points A and F of the region AFED , are defined over the top face of this region. For a reason explained in 
[42], the dual flux (denoted by F(A; F ) or F(F;A)) associated with these two fluxes is defined as the simple 
average of them, i.e., 

F(A; F) = F(F; A) d = [F(AF) + F(FA)] /2 (3.23) 

Similarly, we have 

F(A; F) = F(F; A) = f \P(AF) + F(FA) 1 /2 (3.24) 

The solution-coupling step will be constructed so that, at the top face of any BCE (such as the space-time 
region AFED) sandwiched between the nth and (n — l/2)th time levels, the dual flux evaluated using u r j' 
and ( Ux)(j , (j, n) G fi, is equal to that evaluated using ft" and (h$)”, (j, n) G fi, i.e., 

F(A-F) = F(A-F) (j»Gfi (3.25) 

By substituting Eqs. (3.6), (3.8), (3.15), and (3.17) into Eqs. (3.23) and (3.24), and then imposing the 
condition Eq. (3.25), one has 

( u +) 7 j + ( u -)j+ 1/ 2 = (“+)/ + (w_)" + i/2 0> n ) € ^ (3.26) 

Recall that: (i) (u + )" and (u_)" +1 y 2 , respectively, are the final values of u at the midpoint M + evaluated 
using the solution data stored at the two hosts (points A and F) of the region AFED-, and (ii) (u + )" and 
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(u_)" +1 /2 , respectively, are the intermediate-value counterparts of (u+)" and (u_)" +1 , 2 . As such, Eq. (3.26) 
states that the average of the two final values of u at the midpoint M + evaluated, respectively, using the 
solution data stored at the two hosts of the region AFED, is identical to its intermediate-value counterpart. 
Let 

(«+)" = (1 - V)(u+)] + V («-)"+i /2 (j, n) G 0 (3.27) 

where 77 is a real adjustable parameter. Then (u+)” is an weighted average of (u+)" and (u_)” +1 y 2 - In turn, 
by substituting Eq. (3.27) into Eq. (3.26), one has 

( u -)j+ 1/2 = ? 7(«+)j + (! - »?)(&- )"+i/2 (j» G II (3-28) 

i.e., (u-)] +1/2 is also an weighted average of («+)” and (w_)" +1 , 2 . In fact, Eqs. (3.27) and (3.28) =>• (i) 

(«+)" = ('«+)” and (u_)" +1/2 = («-)7 +1/2 if 17 = 0 (3.29) 

(ii) 

(u+)j = {u-)]+i/2 = \ («+)” + («-)"+i/2 if V — 1/2 (3.30) 

and (iii) 

(«+)" = CMj+i/2 and (m-)"+i / 2 = («+)" if ? 7= 1 (j,n)GO (3.31) 

Because (j, n) G S2 <t=> (j — 1/2, n) G fl, Eq. (3.27) -(3.31) remain valid if each index j in these equations 
is replaced by j — 1/2. In particular, Eq. (3.28) <t=> 

(«-)" = V (m+)"-i/ 2 + (f - »?)(«-)" (/. n) G O (3.32) 


In turn, by substituting Eqs. (3.7), (3.9), (3.12), (3.14), (3.18), and (3.22) into Eqs. (3.27) and (3.32), one 
has 


(u + %)" = (1 - 1l){u + Ux)j + - Ux)j + l/2 

(j, n) G 0 

(3.33) 

{' U - Ui)" = (1 - r])(u - %)" + 77 (u + u $ )"_ 1/2 

(j, n) G 0 

(3.34) 

By adding Eqs. (3.33) and (3.34) together, one has 



«" = (!- + fa/2) (U - «x)"+l/2 + (* + «x)"-l/2 

(j, n) G 0 

(3.35) 


On the other hand, by subtracting Eq. (3.34) from Eq. (3.33), one has 

(«*)" = (1 - V)(ux)j + (V 2 ) (u~ «x)"+i/2 - (u + «s)"_i/2 (j,n)eQ (3.36) 

Eqs. (3.35) and (3.36) form the solution coupling step, i.e., with the aid of these two equations, at each mesh 

point (j, n) G fl, u" and (u $ )" can be determined in terms of the intermediate mesh values at the mesh 
points (j - 1/2, n), (j,n), and (j + 1/2, n). 

As expected from Eq. (3.29), we have 

u] = it] and (u s )" = (u s )" if 77 = 0 (j, n) G O (3.37) 

i.e., solution coupling does not occur when 17 = 0. On the other hand, because of Eq. (3.30), the case with 

17 = 1/2 is referred as the full coupling case. 
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Let the dual-scheme marching step be executed using the a-e scheme, i.e. , Eqs. (3.3) and (3.4) are 
assumed. Then, by substituting Eqs. (3.3) and (3.4) into Eqs. (3.35) and (3.36), one obtains 

u] = | [( 2 - e)u - vux) n r 1/2 

+ ~ — 2 ”^! v ^ u — v ~) u x\ j-i/2 + [(-*- — v ) u ~ (-*- — h ' 2 ) u x\ j + i/2 } Cb n ) e ^ (3.38) 

+ |{[( e + v ) u + ( 2e + v ~ + [( e - v ) u + ( l ’+ v 2 - 2e ) u s]™ + i /2 } 

and 

(«*)" = | D'w + ( 2 - 2e - ^ 2 )u $ ]" _1/2 

+ — 2 ^ { [(e - 1 )u + (2e - 1 + + [(1 - e)« + (2e - 1 - ^)«x] (•?’> n ) e 0 ( 3 - 39 ) 

- |{ [(e + i/) w + (2e + v- i' 2 )ux]J_* /2 - [(e - v)u + (v + i/ 2 - 2e)w $ ]” +1 1/ " j 

Because Eqs. (3.38) and (3.39) are derived assuming that (i) the dual-mesh marching step is carried out using 
the dual a-e scheme, and (ii) the solution-coupling step is formed using an weighted-averaging procedure 
involving the parameter (see Eq. (3.27)), hereafter the full solution-coupling scheme formed by Eqs. (3.38) 
and (3.39) will be referred to as the a-e-r] scheme. 

Alternatively, the dual marching step can be executed using the dual c-r scheme which are defined 
by Eqs. (1.12) and (1.16). For this case, the intermediate values u" and (w s )" should be evaluated using 
Eq. (3.3) and 

(%)" = 2(1 + t) ( ^ ^ + 2v ~ T ) u ^j+i/2 - [u + (1 - 2u - t)u s ]^ZI /2 } (J, n) € ft (3.40) 

By substituting Eqs. (3.3) and (3.40) into Eqs. (3.35) and (3.36), one obtains the c-r-r ] scheme, i.e., 

u i = 2(1 + r ) t (2 + r)u - 2 ^ir 1/2 

+ — 9 + v)u + (1 - + [(1 - v)u - (1 - iy 2 )u s ]^~l^\ 

21 3 ' 3 /1 (.(■") C <» (3.41) 

+ 4(1 + ^ { t( r + v + Tl/ ) u + [2(^ + T ) - (1 + r)v 2 ]u s ] ™_l /2 
+ [(t-v- tv)u + [2(v - t) + (1 + t)v 2 ]u s \ " +1 1/2 | 

and 

(«*)" = ^^ 7 ) K 1 + T ) i/U + I 2 - (! + T ) v ' 2 ] u x] n ~ 1/2 

- 2 h + r) { [“ + (! - 2v - T ) u ■*] - [« - (1 + 2i/ - r)u s ] } 

^ ' (j.n)efi (3.42) 

- 4 ( 1 + r ) {[( T + iy + rt/ ) u + [ 2 ( I/ + ' r ) “ (i + r)^ 2 ]^];- 172 

- [(r - z/ - tv)u + [2(v - r) + (1 + t)u 2 ]u^\ J”* 72 } 
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The von Neumann stability analysis was carried out for both a-e - 77 and c-T-r] scheme. In particular, 
the amplification matrices and the associated amplification factors were obtained for both schemes involving 
arbitrary values of e, r, rj, and the phase angle 9 [72], Because the a-e - 77 and c-r -77 schemes reduce to 
the dual a-e and c-r schemes, respectively when 77 = 0 , as expected, the stability conditions of these two 
schemes also reduce to those of the a-e and c-t schemes, respectively when 77 = 0. For the general case 
involving an arbitrary value of 77 , it turns out that the stability conditions of these schemes generally become 
so complicated that they can only be expressed numerically. As such, the discussion of their stability 
conditions for the general case is beyond the scope of the current paper. 

However, for the full coupling case, i.e., when 77 = 1/2, it can be shown that [72] the a-e-rj scheme is von 
Neumann stable <t=> 

is 4 < e < 1 (the a-e-rj scheme with 77 = 1/2) (3.43) 

Let 

e = /3 M (3> 0 (3.44) 

where (3 > 0 is an adjustable parameter. Then Eq. (3.43) becomes 

v 4 <PW\<l (3-45) 


Excluding the meaningless case v = 0, Eq. (3.45) <t=> (| i ^| 3 — (3) < 0 and \u\ < 1 //?, i.e., 

'1 if 0 = 1 

|^| < minj/3 1 / 3 , 1//3} = <( /3 1/3 < 1 if 0 < (3 < 1 ( 1 / ^ 0) 


(3.46) 


.1/(3 <1 if (3 > 1 


In turn, Eqs. (3.44) and (3.46) => 

e = (3\v\ = 


(^ 0 ) 


(3.47) 


1 if (3 > 1 

/3 4 / 3 <1 if 0 < (3 < 1 

Also, for the case 77 = 1 / 2 , it has been shown that the stability conditions of the c-r - 77 scheme is 

jV) < r (v 2 < 2) (3.48) 


if v 2 < 2 , where 


Note that (i) 


v sU 1 + 2 ) 

1 + 2 s — s- 


(0 < s < 2) 


1 + 2s-s 2 = -(s + \/2-1)(s-\/2-1) >0 if 0<s<2 


(3.49) 


(3.50) 


i.e., the denominator of the expression on the right side of Eq. (3.49) is not zero throughout the domain of 
/(s) and therefore the function is well defined over its domain; and (ii) 


df(s) 2(2 s 2 + s + 1 ) 


>0 if 0 < s < 2 


ds (1 + 2s — s 2 ) 2 
i.e., /(s) is monotonically increasing in its domain and therefore 

0 = /( 0) < /(s) < /( 2) = 8 (0 < s < 2) 


(3.51) 


(3.52) 
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The stability conditions of the c-T-77 scheme for the case 77 = 1/2 and v 2 > 2 can also be cast into an 
analytical form. However, the form is much more complicated and will not be presented here. 

3.2. 2D and 3D extensions 

For a 2D case using triangular meshes, or a 3D case using tetrahedral meshes, there is a natural extension 
of each of the ID a-e -77 and c-r -77 schemes. As an example, consider Fig. 8 where the triangle A BDF and 
its three neighbors A BFH, A BID, and A DJF (which all lie on the x-y plane) are shown. As explained 
in Sec. 2, each of the centroid mesh points G, A, C, and E shown in Fig. 3 is replaced in Fig. 8 by another 
associated interior point which is (i) defined by Eqs. (2.10) and (2.29) with £ > 1; and (ii) denoted by the 
same symbol and also marked by an open circle. The spatial solution point G' (also marked by a cross) 
associated with A BDF again is defined as the centroid of the hexagon ABCDEF . Other spatial solution 
points A' , C", and E' are defined similarly. In addition, points K, L, and M (marked by triangles) depicted 
in Fig. 8, by definition, are the centroids of the quadrilaterals GFAB, GBCD, and GDEF, respectively. 

Let each spatial solution point be assigned an unique identification index (an integer). As an example, 
let points G ' , A', C' , and E' be identified by the indices j, j 1 , j' 2 , and j 3 , respectively. As such, as an 
example, the x- and y- coordinates of point G' will be denoted by Xj and y.j , respectively. Moreover, the 
space-time solution points (G',n), (A’,n), ( C',n ), and ( E',n ), which lie on the nth time level and have 
the points G', A', C' , and E' being their spatial projections, respectively, will be denoted by (j, n) , (ji,n), 
(j 2 ,n), and (j 3 ,n), respectively. Also each space-time solution point (j, n) is associated with [13] (i) one SE 
denoted by SE(j, n); and (ii) three BCEs denoted by CE( fc )(j, n), k = 1,2,3, respectively, and a CCE which 
is the union of CE( fc ) (j, n), k = 1,2,3, and denoted by CE (j,ri). By definition, CE(!)(j, n) is the space-time 
cylinder sandwiched between the nth and (n — l/2)th time levels with the quadrilateral GFAB being its 
spatial projection. CE( 2 )(j, n) and CE( 3 j(j, n) are defined identically except that the quadrilateral GFAB 
are replaced by the quadrilaterals GBCD and GDEF, respectively. It is seen that (i) the solution points 
(j, n) and (ji,n) are the hosts of the space-time region occupied by CE^^j, n), (ii) the solution points ( j , n) 
and (j 2 ,n) are the hosts of the space-time region occupied by CE( 2 )(j, n), and (iii) the solution points (j, n) 
and (j 3 ,n) are the hosts of the space-time region occupied by CE( 3 )(j, n), 

Let the 2D counterpart of Eq. (1.1) be 


du du du 

yrr + a x y — I - a y y~ = 0 

dt dx dy 


(3.53) 


where a x , and a y are constants. Then the 2D counterparts of Eqs. (1.4) and (1.6) are 


h*(x, y, t ; j, n ) = f [a x u* (x, y, t ; j, n),a v u*(x,y,t;j,n),u*(x,y,t; j, n) 


(3.54) 


and 

u*{x, y, t ; j, n) = u] + ( u x )] [(x - Xj) - a x (t - t n )\ + ( u y )" [(y - yj) - a y (t - t n )\ (3.55) 

respectively. Thus there are three independent marching variables, i.e., u", ( u x )j , and (w y )” associated with 
a space-time solution point (j, n) (i.e., point ( G',n )). 

Next let (i) ( K,n ), ( L,n ), and (M, n) denote the points on the nth time level with their spatial pro- 
jections being points K , L , and M, respectively; (ii) the intermediate solution data stored at (j,n), (ji,n), 
(j 2 ,n), and (j 3 ,n) be evaluated from the given solution data at the (n — l/2)th time level using the 2D 
versions of the a-e or c-r scheme [7,13,31]; (iii) [u(K, ?r)]™ and [u(K, n)]^ denote the intermediate values of u 
at ( K,n ) evaluated using the intermediate solution data stored at the hosts (j, n) and (ji,n) of CE^) (j,n), 
respectively; (iv) [ii(L, n)]" and [u(L, n)]™, denote the intermediate values of u at ( L , n) evaluated using the 
intermediate solution data stored at the hosts (j, n) and (. 72 , 77 ) of CE( 2 )(j, n), respectively; (v) [u(M, ?i)]” 
and [ u(M , n)]" 3 denote the intermediate values of u at ( L , n) evaluated using the intermediate solution data 
stored at the hosts (j, n) and (js,n) of CE( 3 )(j, n), respectively; (vi) [u(K, n)]" and [u(K, n)]]? denote the 
final values of u at ( K,n ) evaluated using the final solution data stored at the hosts (j, n) and (ji,n) of 
CE(!)(j, 77 .), respectively; (vii) [u(L,n)]j and [ u ( L ,? 7)]” 2 denote the final values of u at ( L,n ) evaluated using 
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the final solution data stored at the hosts (j, n) and (ji,n) of CE( 2 )(j, n), respectively; and (viii) [u(M, n)]" 
and [u(M, n)]" 3 denote the final values of u at (M, n) evaluated using the final solution data stored at the 
hosts (j,n) and (j3,n) of CE (3)(j, n), respectively. 

Moreover, for each k = 1,2,3, let (i) A k (j,n ) denote the area of the spatial projection of CE( fe )(j, n); 
(ii) F k (j, n) denote the intermediate value of the flux of h* leaving CE( fc )(j, n) through its top face evaluated 
using the intermediate solution data stored at the space-time solution point (j, n ) ; (iii) Fk(jk,n) denote the 
intermediate value of the flux of h* leaving CE( fc )(j, n) through its top face evaluated using the intermediate 
solution data stored at the space-time solution point (iv) Fk(j,n ) denote the final value of the flux 

of h* leaving CE( fc )(j, n) through its top face evaluated using the final solution data stored at the space-time 
solution point ( j,n ); and (v) F k (jk,n ) denote the final value of the flux of h* leaving CE (fc)(j, n) through its 
top face evaluated using the final solution data stored at the space-time solution point ( jk,n ). According 
to the above definitions, A- t (j, n), A 2 (j,n), and A 3 (j,n) are the areas of the quadrilaterals GFAB , GBCD , 
and GDEF, respectively. Because these quadrilaterals are the spatial projections of CE (fc)(j, n), k = 1,2,3, 
respectively, for each k = 1,2,3, Ak(j,n ) is the area of the top face of CE^fj, n). 

With the aid of (i) the above definitions, (ii) the fact that point K depicted in Fig. 8 is the centroid 
of the quadrilateral GFAB (i.e., the spatial projection of CE(!)(j, n)), and (iii) Eq. (3.55) is a first-order 
Taylor’s expansion, by using Eqs. (3.54) and (3.55), one can establish the following relations for several fluxes 
of h* leaving CE^) (j,n) through its top face: 


A (j, n) = [u{K, n)]j Ai ( j , n) (3.56) 

A (ji , n) = [u(K, n)]” Ax ( j , n) (3.57) 

F i (. 3 , n) = [u(K, n)]"Ai (j, n) (3.58) 

and 

F\ (ji , n) = [u(K, n)]” A x ( j , n) (3.59) 

The current solution-coupling step will be constructed assuming the following dual-flux conservation relation: 

Fi ( j , n ) + F\ ( jx , n) = A {j, n) + A (ji ,n) (3.60) 

By substituting Eqs. (3.56)-(3.59) into Eq. (3.60), one obtains 

[u(K, n)}] + [ u(K , n)\l = [ u(K , n)}” + [u(K, n)]£ (3.61) 


Similarly, by considering the fluxes at the top faces of CE( 2 )(ji n ) an( i CE( 3 )(j, n) and imposing the dual-flux 
conservation relations similar to Eq. (3.60), one has 


[u{L, n)]j + [u(L, n)]j 2 = [ u(L , n)]] + [u(L, n)]] 2 (3.62) 

[u{M, n)]j + [ u(M , n)]” = [u(M, n)]] + [u(M, n)]? (3.63) 

We assume the following weighted-averaging relations: 

[u(K, n)]j = (1 - J?) [u(K, n)}] + r] [u(K, n)]” (3.64) 

[u(L, n)]j = (1 -rj) [ u(L , n)]] + 77 [u(L, n)]] 2 (3.65) 

and 

[u(M, n)]j = (1 - 77) [u(M, n)]" + 77 [u(M, n)}] 3 (3.66) 
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where 77 is a real adjustable parameter. In turn, by substituting Eqs. (3.64)-(3.66) into Eqs. (3.61)-(3.63), 
one has 

[u{K, n)]£ = 77 [u(K, n)]” + (1 - rj) [ u(K , n)]£ (3.67) 

n)]" 2 = 77 [u(A, n)]” + (1 - 77) [w(T, n)]? (3.68) 

and 

[u(M, n)]j 3 = 7 ? [u(M, tt)]” + (1 - 77 ) [u(M, n)]” (3.69) 

respectively. 

According to Eqs. (3.64)-(3.66), [u(K,n)]j, [u(L,n)]” and [u(M, ti)]" can be determined in terms of 
[u(K, n)]”, [u(AT, 77.)]” , [u(L,n)]", (fi(L, 71)]”, [w(M, 77)]”, [u(M,n)]" s , and a given value of 77. On the other 
hand, because t = t n at points (7, n), (K,n), (L,n), and (M,ri), Eq. (3.55) along with the definitions of 
[u(K,n)]y, [■ u(L,n )]!?, and [u(M,n)]!>, => 


+ K)"[a;(A:) - ajj] + (u v )][y{K) - y 3 \ = [u(K,n)]] 


(3.70) 


< + K)?[a;(A) - Sj] + (u y )”[7/(L) - %] = [u{L,n)]J 


(3.71) 


and 


< + (ux)j[x(M) - xj] + (u v )][y(M) - yj] = [u(M,ti)]!? 


(3.72) 


Here (i) x(K) and y(K) are the x- and y- coordinates of point K, (ii) x(L) and y(L) are the x- and y- 
coordinates of point L , and (iii) x(M) and y(M) are the x- and y- coordinates of point M. 

Let 

/ 1 x(K) - Xj y{K) - y 3 \ 

/ x(L) - x(K) y{L) - y(K) 

det | 1 x(L) — Xj y(L) — >u I = det | | ^ 0 (3.73) 


1 x{L) - x 3 y(L) - yj 
V 1 x(M) - Xj y(M) - yj J 


x(M) — x(K) y(M) — y(K) 


i.e., points K , L , and M do not lie on a straight line. Then Eqs. (3.70)-(3.72) can be inverted and, as a 
result, u", (u x )", and (■ u y ) r j can be explicitly determined in terms of [u(K,n)]j, [u(L,r i)]J, and [u(M, n)]” 
if the spatial coordinates of points (j, n) , AT, L, and M are given. As such, given all the spatial geometric 
data, the final solution values at any (j. n) can be determined in terms of the given solution values at the 
(n — 1/2) time level by using the 2D solution-coupling procedure describe above. 

In a CESE 3D scheme using a spatial mesh constructed from tetrahedra, generally a solution point is 
associated with four independent mesh variables, four BCEs and one CCE. Also each BCE has two hosts. As 
such, construction of the 3D solution-coupling procedure is a straightforward extension of its 2D counterpart. 

Finally, it should be emphasized that the use of a solution-coupling procedure generally 
will introduce numerical dissipation and, as a result, degrade accuracy. As such it is strongly 
recommended that, as long as it serves the purpose of preventing solution decoupling, this 
procedure should be applied as sparingly as possible, i.e., it should be used with a small value 
of y and/or applied once after a preset number of marching steps. 


4. Conclusions and discussions 


As a preliminary for the later development, a review of basic CESE ideas and schemes along with a 
discussion of several issues related to the application of the CESE method involving unstructured triangu- 
lar/tetrahedral meshes was presented in Sec. 1. Among the issues discussed are: (A) geometric difficulties 
associated with a boundary which has high curvature or sharp corner and is also surrounded by triangu- 
lar/tetrahedral meshes of extremely high aspect ratio, and (B) possible occurrences of solution decoupling 
in CESE simulations involving unstructured triangular/tetrahedral meshes. 

To tackle the issue (A), Sec. 2 began by identifying the root cause of the geometry difficulty for the 
triangular-mesh case, i.e., the hexagon ABC DEF depicted in Fig. 4 may become highly concave if (i) the 
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mesh points A, C, E, and G are designated, respectively, as the centroids of the associated triangles, and 
(ii) two neighboring triangles with high aspect ratio also have sharply different orientations. A rigorous 
mathematical analysis was then presented to validate the observation. According to this analysis, the 
centroid and in-center of a triangle are but two special cases with £ = 0 and Q = 1, respectively, of the 
family of the interior points defined by Eqs. (2.10), (2.29), and (2.32). By using the results presented in 
Eqs. (2.39)-(2.41), it becomes clear that, by replacing the centroid with another interior point with £ > 1, 
the hexagon ABC DEF referred to above will become convex or much less concave. This conclusion has been 
verified by numerous numerical experiments. In fact, the geometric difficulty can be overcome by simply 
designating points A, C, E, and G (see Fig. 4), respectively, as the in-centers (£ = 1) of the associated 
triangles. Moreover, it was shown in Sec. 2 that the tetrahedral-mesh case can be dealt with in a completely 
parallel manner. 

The new solution-coupling procedures described in Sec. 3 have become an integral part of the CESE 
code development. In fact, the nuisance of solution decoupling is avoided completely with the use of these 
new procedures. 
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S(V) of an arbitrary space-time volume V. 
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2(a) — The space-time mesh. 




(j— 1/2,n— 1/2) ° 0+1/2, n-1/2) 
2(e) CE(j,n) 

Figure 2. — The SEs and CEs. 
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Figure 3. — A spatial domain formed by triangles. 
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(xi, yi) 



Figure 5. — The perpendicular projections of paint Po on the 
three sides (or their extensions) of a triangle AP-| , P 2 , P 3 . 
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(j— 1/2, n— 1/2) 


(j, n— 1/2) 


0+1/2, n— 1/2) 



Figure 7. — Construction of the a-e-rj and c-x-r| schemes. 
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